---
title: "saliencyVAS"
output: html_document
date: "2024-09-24"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(base)       # The R Base Package
library(rstudioapi) #to set directory to current file path
library(matrixStats)# Functions that Apply to Rows and Columns of Matrices (and to Vectors) 
library(sjmisc)     # Data and Variable Transformation Functions 
library(stringr)    # Simple, Consistent Wrappers for Common String Operations (filtering of partial name)
library(purrr)      # Functional Programming Tools

library(ggplot2)    # Create Elegant Data Visualizations Using the Grammar of Graphics
library(ggpubr)     # 'ggplot2' Based Publication Ready Plots
library(ggsignif)   # Significance Brackets for 'ggplot2'
library(ggsci)      # Scientific Journal and Sci-Fi Themed Color Palettes for 'ggplot2'
library(ggrain)     # A Rainclouds Geom for 'ggplot2'
library(plotly)     # smoothens the lines in plots, especially the horizontal ones
library(RColorBrewer) # color palette 
cbPalette <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")# color-blind friendly palette

library(psych)      # Procedures for Psychological, Psychometric, and Personality Research
library(irr)        # Various Coefficients of Interrater Reliability and Agreement
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests 
library(broom)      # Convert Statistical Objects into Tidy Tibbles 
library(stats)      # The R Stats Package
library(lme4)       # Linear Mixed-Effects Models using 'Eigen' and S4 
library(lmerTest)   #get p-values for LMM's
library(emmeans)    # Estimated Marginal Means, aka Least-Squares Means
library(car)        # Companion to Applied Regression
library(pbkrtest)   # Parametric Bootstrap, Kenward-Roger and Satterthwaite Based Methods for Test in Mixed Models
library(effectsize) # Indices of Effect Size

library(sjlabelled) # to generate LMM output table 
library(sjPlot)     # table functions
```

### set working directory
```{r}
current_path = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(current_path ))
print( getwd() )
```

### load data
```{r, include=FALSE}
highOddball <-read.delim("VAS_rating_peaksHighO.csv", header=TRUE, sep=";")
lowOddball <-read.delim("VAS_rating_peaksLowO.csv", header=TRUE, sep=";")
HO_base <-read.delim("VAS_rating_peakHighBase.csv", header=TRUE, sep=";")
LO_base <-read.delim("VAS_rating_peakLowBase.csv", header=TRUE, sep=";")
```

### sort table
```{r}
highOddball <-convert_as_factor(highOddball, vars=  c("condition", "stimulation","subject"))
lowOddball <-convert_as_factor(lowOddball, vars=  c("condition", "stimulation","subject"))
HO_base <-convert_as_factor(HO_base, vars=  c("condition", "stimulation","subject"))
LO_base <-convert_as_factor(LO_base, vars=  c("condition", "stimulation","subject"))
```

```{r}
all_peaks <- rbind(highOddball, lowOddball, HO_base, LO_base)
peaks_highO <- rbind(highOddball, HO_base)
peaks_lowO <- rbind(lowOddball, LO_base)
```

### testing LMM assumptions
```{r}
shapiro.test(peaks_lowO$rating)
levene_test(rating~stimulation, data=peaks_lowO)

shapiro.test(peaks_highO$rating)
levene_test(rating~stimulation, data=peaks_highO)

check_data <- na.omit(peaks_highO)
lmm <- lmer(rating ~ stimulation +(1|subject), data=check_data)
qqnorm(resid(lmm))
qqline(resid(lmm))

check_data$cook <- cooks.distance(lmm)
plot(check_data$cook, check_data$subject)

check_data$amp_log <- log(check_data$rating)
lmm_log <- lmer(amp_log ~ stimulation+ (1|subject), data=check_data)
qqnorm(resid(lmm_log))
qqline(resid(lmm_log))

check_data$cook <- cooks.distance(lmm_log)
plot(check_data$cook, check_data$subject)

check_data$amplitude <- replace(check_data$amplitude, c(check_data$cook > 1),NA) #replaces amplitude with NA for all rows in which cook's distance is larger than 1
check_data$amp_log <- replace(check_data$amp_log, c(check_data$cook > 1),NA)

peaks_lowO <- check_data
```

### oddball vs baseline for each condition (LMM)
```{r}
ggboxplot(all_peaks, x="stimulation", y="rating", facet.by = "condition", group="subject", add=c("mean_se"), title= "VAS", xlab= FALSE, ylab="VAS [0-10]")  + geom_exec(geom_line, all_peaks, x="stimulation", y="rating", group="subject", color="grey" ) + stat_compare_means(method="anova")

## stats
LMM_rating <- lmer(rating ~ stimulation+ (1|subject), data=peaks_highO, REML=TRUE)
summary(LMM_rating)
if(requireNamespace("pbkrtest", quietly = TRUE))
anova(LMM_rating, ddf="Kenward-Roger", type =3)

LMM_rating <- lmer(rating ~ stimulation+ (1|subject), data=peaks_lowO, REML=TRUE)
summary(LMM_rating)
if(requireNamespace("pbkrtest", quietly = TRUE))
anova(LMM_rating, ddf="Kenward-Roger", type =3)
```

### relative amplitude difference between oddballs
```{r}
highOddball$relativeVASHigh<- highOddball$rating-HO_base$rating
lowOddball$relativeVASLow<-lowOddball$rating-LO_base$rating
t.test(highOddball$relativeVAS,lowOddball$relativeVAS, paired=TRUE )
wilcox.test(highOddball$relativeVAS,lowOddball$relativeVAS, paired=TRUE )

relativeVASHigh<- as.data.frame(highOddball$rating-HO_base$rating)
colnames(relativeVASHigh) <- "relative rating high oddball"
relativeVASLow<- as.data.frame(lowOddball$rating-LO_base$rating)
colnames(relativeVASLow) <- "relative rating low oddball"

rel_rat <- stack(cbind(relativeVASHigh,relativeVASLow))
peaks_highO_rel <- cbind(peaks_highO, rel_rat)

ggerrorplot(rel_rat, y="values", x="ind", add="jitter")

relRatPlot<-   ggerrorplot(peaks_highO_rel, x="ind", y="values", desc.stat= "mean_ci", combine = TRUE, size= 1, 
          add=c( "jitter"),add.params = list(color="black", alpha=0.5), title= FALSE, xlab= FALSE, ylab="∆ VAS ratings", color="ind", palette = c("Red", "Blue"), legend="RIGHT",
          panel.labs.background=list(color="white", fill="white"), panel.labs.font=list(face="bold", size=11))+
  theme(panel.border = element_blank(), axis.line = element_line()) + 
    scale_x_discrete(labels=c('high oddball', 'low oddball'),expand=expansion(mult=c(0,2)))+
  stat_summary(data= peaks_highO_rel,  fun = "mean", geom = "line", aes(group = subject),color="grey", alpha=0.5) +   stat_compare_means(method="t.test", paired=TRUE, comparisons = list(c("relative rating high oddball", "relative rating low oddball"))  ,label = "p.signif", tip.length = 0.02)

```



```{r}
ggplot(all_peaks, aes(stimulation, rating, fill=condition))+
  geom_rain(alpha=.5,rain.side = "f2x2", cov="condition",
            boxplot.args = list(color="black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15), 
            violin.args.pos = list( side = "r",
            width = 0.7, position = position_nudge(x = 0.3)))+
  scale_colour_manual(values=cbPalette)+
  scale_fill_manual(values=cbPalette)+
  stat_summary(data= peaks_highO,  fun = "mean", geom = "line", aes(group = subject, color=condition),alpha=0.5)+ #alpha=x to make lines transparent
  stat_summary(data= peaks_lowO, fun = "mean", geom = "line", aes(group = subject, color=condition),alpha=0.5)+
  stat_summary(data= peaks_highO, fun = "mean", geom = "line", aes(group = condition, color=condition),size=1.25)+
  stat_summary(data= peaks_lowO, fun = "mean", geom = "line", aes(group = condition, color=condition), size=1.25)+
  scale_y_continuous(breaks = c(2, 4, 6,8,10))+
  theme_classic() +
  labs(x="", y="VAS rating [0-10]")+
  geom_bracket(data=peaks_highO, xmin = c("oddball"), xmax = c("baseline"), y.position = c(10.5), label = c("***"),  tip.length = 0.00,aes(color=condition), vjust = 0.5, bracket.shorten = 1.8, size = 1)+
    geom_bracket(data=peaks_lowO, xmin = c("oddball"), xmax = c("baseline"), y.position = c(11), label = c("ns"),  tip.length = 0.00,aes(color=condition),  vjust = 0.1, bracket.shorten = 1.5,  size = 1)+
theme(legend.title = element_text(size=14,face = "bold"), legend.text = element_text(size=13), axis.text.x = element_text(size=14, vjust=0.5), axis.text.y = element_text(hjust = 0.50,size=14), axis.title.x = element_text( face = "bold",size=14), axis.title.y = element_text(face = "bold", size=14))

```
